IGCCodonSimulationSummary

This script is used for summarizing IGC Simulation result (with tract length)

Read in data first

rm(list=ls())  # clean up workspace
IGC.path <- "/Users/xji3/FromCluster_IGCCodonSimulation01262016/"
paml.path <- "/Users/xji3/GitFolders/IGCCodonSimulation/"
IGC.geo.list <- c(3.0, 10.0, 50.0, 100.0, 500.0)
#IGC.geo.list <- c(1.0)
num.sim = 100
paralog = "YDR418W_YEL054C"


# Now read in PAML analysis of these simulated data set
data.path <- paste(paralog, "",sep = "/")
for (IGC.geo in IGC.geo.list){
  summary_mat <- NULL
  IGC.geo.path <- paste("IGCgeo_", toString(IGC.geo), ".0/", sep = "")
  file.name <- paste("geo", paste(toString(IGC.geo), ".0", sep = ""), "paml", "unrooted", "summary.txt", sep = "_")
  for (sim.num in 0:(num.sim - 1)){
    summary_file <- paste(paml.path, file.name, sep = "")
    if (file.exists(summary_file)){
      all <- readLines(summary_file, n = -1)
      col.names <- strsplit(all[1], ' ')[[1]][-1]
      row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
      summary_mat <- as.matrix(read.table(summary_file, 
                                          row.names = row.names, 
                                          col.names = col.names))
      
      }
    }
  assign(paste("PAML", paste(toString(IGC.geo), ".0", sep = ""), "summary", sep = "_"), summary_mat)
  }

#Read.summary <- function(IGC.path, IGC.x3.path, IGC.geo.list, num.sim = 100, paralog = "YDR418W_YEL054C"){
## Now read simulation result

sim.path <- paste("SimulationSummary", paralog, "",sep = "/")
for (IGC.geo in IGC.geo.list){
  summary_mat <- NULL
  IGC.geo.path <- paste("IGCgeo_", toString(IGC.geo), ".0/", sep = "")
  file.prefix <- paste(paralog, "MG94_geo", paste(toString(IGC.geo), ".0", sep = ""), "Sim", sep = "_")
  file.suffix <- "summary.txt"
  for (sim.num in 0:(num.sim - 1)){
    file.name <- paste(file.prefix, toString(sim.num), file.suffix, sep = "_")
    summary_file <- paste(IGC.path, sim.path, IGC.geo.path, file.name, sep = "")
    if (file.exists(summary_file)){
      all <- readLines(summary_file, n = -1)
      col.names <- paste("geo", paste(toString(IGC.geo), ".0", sep = ""), "sim", toString(sim.num), sep = "_")
      row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
      summary_mat <- cbind(summary_mat, as.matrix(read.table(summary_file, 
                                                             row.names = row.names, 
                                                             col.names = col.names)))        
      }
    }
  assign(paste("geo", paste(toString(IGC.geo), ".0", sep = ""), "summary", sep = "_"), summary_mat)
  }


# Now read in 'real'  % changes due to IGC in the simulation
data.path <- paste(paralog, "",sep = "/")
for (IGC.geo in IGC.geo.list){
  summary_mat <- NULL
  IGC.geo.path <- paste("IGCgeo_", toString(IGC.geo), ".0/", sep = "")
  file.prefix <- paste(paralog, "MG94_geo", paste(toString(IGC.geo), ".0", sep = ""), "Sim", sep = "_")
  file.suffix <- "short.log"
  for (sim.num in 0:(num.sim - 1)){
    file.name <- paste(file.prefix, toString(sim.num), file.suffix, sep = "_")
    file.name <- paste("sim", paste(toString(sim.num), "/", file.prefix, sep = ""),
                       toString(sim.num), file.suffix, sep = "_")
    summary_file <- paste(IGC.path, data.path, IGC.geo.path, file.name, sep = "")
    if (file.exists(summary_file)){
      all <- readLines(summary_file, n = -1)
      col.names <- paste("geo", paste(toString(IGC.geo), ".0", sep = ""), "sim", toString(sim.num), sep = "_")
      row.names <- strsplit(all[length(all)], ' ')[[1]][-1]
      summary_mat <- cbind(summary_mat, as.matrix(read.table(summary_file, 
                                                             row.names = row.names, 
                                                             col.names = col.names)))        
      }
    }
  assign(paste("pIGC", paste(toString(IGC.geo), ".0", sep = ""), "summary", sep = "_"), summary_mat)
  }

Now analyze the simulation results

parameters used in generating datasets are:

pi_a = 0.2983654, pi_c = 0.2095729, pi_g = 0.1871536, pi_t = 0.3049081

kappa = 8.4043336

omega = 1.0

tau = 1.409408 (0.42)

Tau estimate summary

# geo_3.0
summary(geo_3.0_summary["tau", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.815   1.200   1.290   1.360   1.550   1.950

sd(geo_3.0_summary["tau", ])

## [1] 0.2547

# geo_10.0
summary(geo_10.0_summary["tau", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.539   1.200   1.450   1.470   1.670   2.650

sd(geo_10.0_summary["tau", ])

## [1] 0.4084

# geo_50.0
summary(geo_50.0_summary["tau", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.261   1.020   1.380   1.490   1.930   3.890

sd(geo_50.0_summary["tau", ])

## [1] 0.67

# geo_100.0
summary(geo_100.0_summary["tau", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.114   1.090   1.600   1.620   1.940   4.500

sd(geo_100.0_summary["tau", ])

## [1] 0.779

# geo_500.0
summary(geo_500.0_summary["tau", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.78    1.78    1.98    2.91    7.14

sd(geo_500.0_summary["tau", ])

## [1] 1.52

sd.list <- c(sd(geo_3.0_summary["tau", ]), 
             sd(geo_10.0_summary["tau", ]),sd(geo_50.0_summary["tau", ]),
             sd(geo_100.0_summary["tau", ]), sd(geo_500.0_summary["tau", ]))
plot(IGC.geo.list, sd.list)
abline(h = 0.42)

mean.list <- c(mean(geo_3.0_summary["tau", ]),
               mean(geo_10.0_summary["tau", ]), mean(geo_50.0_summary["tau", ]), 
               mean(geo_100.0_summary["tau", ]), mean(geo_500.0_summary["tau", ]))
plot(IGC.geo.list, mean.list)
abline(h = 1.409408)

Now summarized % changes due to IGC in simulation estimates

It was estimated 0.2131 (0.02884)

percent.IGC.geo.3.0 <- colSums(geo_3.0_summary[34:45, ] + geo_3.0_summary[46:57, ]) / (colSums(geo_3.0_summary[34:45, ] + geo_3.0_summary[46:57, ] + geo_3.0_summary[58:69, ]))
summary(percent.IGC.geo.3.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.151   0.194   0.209   0.210   0.227   0.265

sd(percent.IGC.geo.3.0)

## [1] 0.02406

percent.IGC.geo.10.0 <- colSums(geo_10.0_summary[34:45, ] + geo_10.0_summary[46:57, ]) / (colSums(geo_10.0_summary[34:45, ] + geo_10.0_summary[46:57, ] + geo_10.0_summary[58:69, ]))
summary(percent.IGC.geo.10.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.133   0.193   0.215   0.215   0.243   0.287

sd(percent.IGC.geo.10.0)

## [1] 0.03418

percent.IGC.geo.50.0 <- colSums(geo_50.0_summary[34:45, ] + geo_50.0_summary[46:57, ]) / (colSums(geo_50.0_summary[34:45, ] + geo_50.0_summary[46:57, ] + geo_50.0_summary[58:69, ]))
summary(percent.IGC.geo.50.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0678  0.1780  0.2140  0.2130  0.2580  0.3750

sd(percent.IGC.geo.50.0)

## [1] 0.05983

percent.IGC.geo.100.0 <- colSums(geo_100.0_summary[34:45, ] + geo_100.0_summary[46:57, ]) / (colSums(geo_100.0_summary[34:45, ] + geo_100.0_summary[46:57, ] + geo_100.0_summary[58:69, ]))
summary(percent.IGC.geo.100.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0286  0.1690  0.2320  0.2230  0.2670  0.3880

sd(percent.IGC.geo.100.0)

## [1] 0.06909

percent.IGC.geo.500.0 <- colSums(geo_500.0_summary[34:45, ] + geo_500.0_summary[46:57, ]) / (colSums(geo_500.0_summary[34:45, ] + geo_500.0_summary[46:57, ] + geo_500.0_summary[58:69, ]))
summary(percent.IGC.geo.500.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.137   0.232   0.226   0.326   0.488

sd(percent.IGC.geo.500.0)

## [1] 0.1272

# plot mean % changes due to IGC v.s. IGC_geo
mean.percent.IGC.list <- c( mean(percent.IGC.geo.3.0), 
                            mean(percent.IGC.geo.10.0), mean(percent.IGC.geo.50.0), 
                            mean(percent.IGC.geo.100.0), mean(percent.IGC.geo.500.0))
plot(IGC.geo.list, mean.percent.IGC.list)
abline(h = 0.2131)

sd.percent.IGC.list <- c(sd(percent.IGC.geo.3.0), 
                         sd(percent.IGC.geo.10.0), sd(percent.IGC.geo.50.0), 
                         sd(percent.IGC.geo.100.0), sd(percent.IGC.geo.500.0))

plot(IGC.geo.list, sd.percent.IGC.list)
abline(h = 0.02884)

Now get histograms of tau estimate

# geo_3.0
hist(geo_3.0_summary["tau",])

# geo_10.0
hist(geo_10.0_summary["tau",])

# geo_50.0
hist(geo_50.0_summary["tau",])

# geo_100.0
hist(geo_100.0_summary["tau",])

# geo_500.0
hist(geo_500.0_summary["tau",])

Now get histograms of % change from IGC from estimate

# geo_3.0
hist(percent.IGC.geo.3.0)

# geo_10.0
hist(percent.IGC.geo.10.0)

# geo_50.0
hist(percent.IGC.geo.50.0)

# geo_100.0
hist(percent.IGC.geo.100.0)

# geo_500.0
hist(percent.IGC.geo.500.0)

Now summarize real % change due to IGC in simulation

It was estimated 0.2131 (0.02884) in real data set

pIGC.real.geo.3.0 <- pIGC_3.0_summary[1, ]
summary(pIGC.real.geo.3.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.168   0.200   0.213   0.215   0.228   0.262

sd(pIGC.real.geo.3.0)

## [1] 0.01977

pIGC.real.geo.10.0 <- pIGC_10.0_summary[1, ]
summary(pIGC.real.geo.10.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.139   0.198   0.214   0.215   0.238   0.279

sd(pIGC.real.geo.10.0)

## [1] 0.02781

pIGC.real.geo.50.0 <- pIGC_50.0_summary[1, ]
summary(pIGC.real.geo.50.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0698  0.1790  0.2140  0.2120  0.2520  0.3150

sd(pIGC.real.geo.50.0)

## [1] 0.0522

pIGC.real.geo.100.0 <- pIGC_100.0_summary[1, ]
summary(pIGC.real.geo.100.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0337  0.1790  0.2170  0.2170  0.2530  0.3580

sd(pIGC.real.geo.100.0)

## [1] 0.06202

pIGC.real.geo.500.0 <- pIGC_500.0_summary[1, ]
summary(pIGC.real.geo.500.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.146   0.214   0.210   0.273   0.436

sd(pIGC.real.geo.500.0)

## [1] 0.1049

# plot mean % changes due to IGC v.s. IGC_geo
mean.pIGC.real.list <- c( mean(pIGC.real.geo.3.0), 
                          mean(pIGC.real.geo.10.0), mean(pIGC.real.geo.50.0), 
                          mean(pIGC.real.geo.100.0), mean(pIGC.real.geo.500.0))
plot(IGC.geo.list, mean.pIGC.real.list)
abline(h = 0.2131)

sd.pIGC.real.list <- c(sd(pIGC.real.geo.3.0), 
                       sd(pIGC.real.geo.10.0), sd(pIGC.real.geo.50.0), 
                       sd(pIGC.real.geo.100.0), sd(pIGC.real.geo.500.0))

plot(IGC.geo.list, sd.pIGC.real.list)
abline(h = 0.02884)

Now get histograms of real % change due to IGC from simulation

# geo_3.0
hist(pIGC.real.geo.3.0)

# geo_10.0
hist(pIGC.real.geo.10.0)

# geo_50.0
hist(pIGC.real.geo.50.0)

# geo_100.0
hist(pIGC.real.geo.100.0)

# geo_500.0
hist(pIGC.real.geo.500.0)

Now summarize difference between real % change due to IGC with estimated %

# geo_3.0
diff.pIGC.geo.3.0 <- percent.IGC.geo.3.0 - pIGC.real.geo.3.0
summary(diff.pIGC.geo.3.0)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.03820 -0.01680 -0.00743 -0.00451  0.00837  0.03130

sd(diff.pIGC.geo.3.0)

## [1] 0.01688

# geo_10.0
diff.pIGC.geo.10.0 <- percent.IGC.geo.10.0 - pIGC.real.geo.10.0
summary(diff.pIGC.geo.10.0)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.05680 -0.01590  0.00092 -0.00011  0.01410  0.05400

sd(diff.pIGC.geo.10.0)

## [1] 0.02105

# geo_50.0
diff.pIGC.geo.50.0 <- percent.IGC.geo.50.0 - pIGC.real.geo.50.0
summary(diff.pIGC.geo.50.0)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.05880 -0.01320 -0.00050  0.00039  0.01450  0.07750

sd(diff.pIGC.geo.50.0)

## [1] 0.02502

# geo_100.0
diff.pIGC.geo.100.0 <- percent.IGC.geo.100.0 - pIGC.real.geo.100.0
summary(diff.pIGC.geo.100.0)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.04790 -0.01170  0.00460  0.00546  0.01700  0.08940

sd(diff.pIGC.geo.100.0)

## [1] 0.02377

# geo_500.0
diff.pIGC.geo.500.0 <- percent.IGC.geo.500.0 - pIGC.real.geo.500.0
summary(diff.pIGC.geo.500.0)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.06150 -0.01080  0.00762  0.01600  0.04040  0.11700

sd(diff.pIGC.geo.500.0)

## [1] 0.03958

# plot mean % changes due to IGC v.s. IGC_geo
mean.diff.pIGC.list <- c(mean(diff.pIGC.geo.3.0), 
                         mean(diff.pIGC.geo.10.0), mean(diff.pIGC.geo.50.0), 
                         mean(diff.pIGC.geo.100.0), mean(diff.pIGC.geo.500.0))
plot(IGC.geo.list, mean.diff.pIGC.list)
abline(h = 0.0)

sd.diff.pIGC.list <- c(sd(diff.pIGC.geo.3.0), 
                       sd(diff.pIGC.geo.10.0), sd(diff.pIGC.geo.50.0), 
                       sd(diff.pIGC.geo.100.0), sd(diff.pIGC.geo.500.0))

plot(IGC.geo.list, sd.diff.pIGC.list)
abline(h = 0.02884)

Now get histograms of difference of % IGC from ‘real’ and estimate

# geo_3.0
hist(diff.pIGC.geo.3.0)

# geo_10.0
hist(diff.pIGC.geo.10.0)

# geo_50.0
hist(diff.pIGC.geo.50.0)

# geo_100.0
hist(diff.pIGC.geo.100.0)

# geo_500.0
hist(diff.pIGC.geo.500.0)

update Jan 29th, 2016

Now added PAML estimate of all simulated datasets. Start to show more comparisons.

Omega estimate from PAML

# geo_3.0
summary(PAML_3.0_summary["omega", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.677   0.885   1.020   1.010   1.110   1.290

sd(PAML_3.0_summary["omega", ])

## [1] 0.138

# geo_10.0
summary(PAML_10.0_summary["omega", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.605   0.914   0.971   0.996   1.100   1.480

sd(PAML_10.0_summary["omega", ])

## [1] 0.1516

# geo_50.0
summary(PAML_50.0_summary["omega", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.725   0.933   1.020   1.020   1.120   1.490

sd(PAML_50.0_summary["omega", ])

## [1] 0.1596

# geo_100.0
summary(PAML_100.0_summary["omega", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.688   0.904   0.987   0.996   1.060   1.500

sd(PAML_100.0_summary["omega", ])

## [1] 0.1447

# geo_500.0
summary(PAML_500.0_summary["omega", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.606   0.896   1.010   1.020   1.120   1.510

sd(PAML_500.0_summary["omega", ])

## [1] 0.1677

PAML.omega.mean <- c(mean(PAML_3.0_summary["omega", ]), mean(PAML_10.0_summary["omega", ]),
                     mean(PAML_50.0_summary["omega", ]), mean(PAML_100.0_summary["omega", ]), 
                     mean(PAML_500.0_summary["omega", ]))
PAML.omega.sd <- c(sd(PAML_3.0_summary["omega", ]), sd(PAML_10.0_summary["omega", ]),
                   sd(PAML_50.0_summary["omega", ]), sd(PAML_100.0_summary["omega", ]), 
                   sd(PAML_500.0_summary["omega", ]))

plot(IGC.geo.list, PAML.omega.mean)
abline(h = 1.0)

plot(IGC.geo.list, PAML.omega.sd)

Omega estimate from IGC + MG94

# geo_3.0
summary(geo_3.0_summary["omega", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.701   0.921   1.000   1.010   1.090   1.340

sd(geo_3.0_summary["omega", ])

## [1] 0.1234

# geo_10.0
summary(geo_10.0_summary["omega", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.730   0.917   0.984   1.000   1.080   1.430

sd(geo_10.0_summary["omega", ])

## [1] 0.1407

# geo_50.0
summary(geo_50.0_summary["omega", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.748   0.916   1.010   1.020   1.110   1.380

sd(geo_50.0_summary["omega", ])

## [1] 0.1436

# geo_100.0
summary(geo_100.0_summary["omega", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.754   0.914   1.010   1.000   1.060   1.420

sd(geo_100.0_summary["omega", ])

## [1] 0.1314

# geo_500.0
summary(geo_500.0_summary["omega", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.675   0.922   1.010   1.020   1.110   1.460

sd(geo_500.0_summary["omega", ])

## [1] 0.1481

geo.omega.mean <- c(mean(geo_3.0_summary["omega", ]), mean(geo_10.0_summary["omega", ]),
                    mean(geo_50.0_summary["omega", ]), mean(geo_100.0_summary["omega", ]), 
                    mean(geo_500.0_summary["omega", ]))
geo.omega.sd <- c(sd(geo_3.0_summary["omega", ]), sd(geo_10.0_summary["omega", ]),
                  sd(geo_50.0_summary["omega", ]), sd(geo_100.0_summary["omega", ]), 
                  sd(geo_500.0_summary["omega", ]))

plot(IGC.geo.list, geo.omega.mean)
abline(h = 1.0)

plot(IGC.geo.list, geo.omega.sd)

kappa estimate from PAML

# geo_3.0
summary(PAML_3.0_summary["kappa", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.63    7.42    8.44    8.49    9.44   12.80

sd(PAML_3.0_summary["kappa", ])

## [1] 1.373

# geo_10.0
summary(PAML_10.0_summary["kappa", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    6.10    7.58    8.11    8.36    9.03   14.30

sd(PAML_10.0_summary["kappa", ])

## [1] 1.39

# geo_50.0
summary(PAML_50.0_summary["kappa", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.85    7.41    8.26    8.55    9.28   13.00

sd(PAML_50.0_summary["kappa", ])

## [1] 1.556

# geo_100.0
summary(PAML_100.0_summary["kappa", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.29    7.33    8.30    8.49    9.53   13.00

sd(PAML_100.0_summary["kappa", ])

## [1] 1.572

# geo_500.0
summary(PAML_500.0_summary["kappa", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.50    7.56    8.44    8.67    9.53   16.00

sd(PAML_500.0_summary["kappa", ])

## [1] 1.703

PAML.kappa.mean <- c(mean(PAML_3.0_summary["kappa", ]), mean(PAML_10.0_summary["kappa", ]),
                     mean(PAML_50.0_summary["kappa", ]), mean(PAML_100.0_summary["kappa", ]), 
                     mean(PAML_500.0_summary["kappa", ]))
PAML.kappa.sd <- c(sd(PAML_3.0_summary["kappa", ]), sd(PAML_10.0_summary["kappa", ]),
                   sd(PAML_50.0_summary["kappa", ]), sd(PAML_100.0_summary["kappa", ]), 
                   sd(PAML_500.0_summary["kappa", ]))

plot(IGC.geo.list, PAML.kappa.mean)
abline(h = 8.4043336)

plot(IGC.geo.list, PAML.kappa.sd)

kappa estimate from IGC + MG94

# geo_3.0
summary(geo_3.0_summary["kappa", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    6.01    7.54    8.18    8.50    9.20   13.50

sd(geo_3.0_summary["kappa", ])

## [1] 1.321

# geo_10.0
summary(geo_10.0_summary["kappa", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.61    7.59    8.21    8.43    8.99   13.10

sd(geo_10.0_summary["kappa", ])

## [1] 1.308

# geo_50.0
summary(geo_50.0_summary["kappa", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    6.25    7.70    8.43    8.65    9.46   13.40

sd(geo_50.0_summary["kappa", ])

## [1] 1.44

# geo_100.0
summary(geo_100.0_summary["kappa", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.42    7.51    8.43    8.59    9.50   13.00

sd(geo_100.0_summary["kappa", ])

## [1] 1.556

# geo_500.0
summary(geo_500.0_summary["kappa", ])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.71    7.65    8.41    8.65    9.47   14.80

sd(geo_500.0_summary["kappa", ])

## [1] 1.613

geo.kappa.mean <- c(mean(geo_3.0_summary["kappa", ]), mean(geo_10.0_summary["kappa", ]),
                    mean(geo_50.0_summary["kappa", ]), mean(geo_100.0_summary["kappa", ]), 
                    mean(geo_500.0_summary["kappa", ]))
geo.kappa.sd <- c(sd(geo_3.0_summary["kappa", ]), sd(geo_10.0_summary["kappa", ]),
                  sd(geo_50.0_summary["kappa", ]), sd(geo_100.0_summary["kappa", ]), 
                  sd(geo_500.0_summary["kappa", ]))

plot(IGC.geo.list, geo.kappa.mean)
abline(h = 8.4043336)

plot(IGC.geo.list, geo.kappa.sd)

Now start to look at branch length estimates

Branch lengths used in simulation:

Branch	blen
N0_N1:	0.0197240946542
N0_kluyveri	0.215682181791
N1_N2	0.20925129872
N1_castellii	0.171684721483
N2_N3	0.0257112589202
N2_bayanus	0.0266075664688
N3_N4	0.0321083243449
N3_kudriavzevii	0.0853588718458
N4_N5	0.024947887926
N4_mikatae	0.0566627496729
N5_cerevisiae	0.0581451177847
N5_paradoxus	0.0218788166581

# N0_N1
IGC.N0.N1.mean.list <- c(mean(geo_3.0_summary["(N0,N1)", ]), mean(geo_10.0_summary["(N0,N1)", ]),
                         mean(geo_50.0_summary["(N0,N1)", ]), mean(geo_100.0_summary["(N0,N1)", ]),
                         mean(geo_500.0_summary["(N0,N1)", ]))
IGC.N0.N1.sd.list <- c(sd(geo_3.0_summary["(N0,N1)", ]), sd(geo_10.0_summary["(N0,N1)", ]),
                       sd(geo_50.0_summary["(N0,N1)", ]), sd(geo_100.0_summary["(N0,N1)", ]),
                       sd(geo_500.0_summary["(N0,N1)", ]))

PAML.N0.N1.paralog1.mean.list <- c(mean(PAML_3.0_summary["N0_N1", ]), 
                                   mean(PAML_10.0_summary["N0_N1", ]),
                                   mean(PAML_50.0_summary["N0_N1", ]),
                                   mean(PAML_100.0_summary["N0_N1", ]),
                                   mean(PAML_500.0_summary["N0_N1", ]))
PAML.N0.N1.paralog1.sd.list <- c(sd(PAML_3.0_summary["N0_N1", ]), 
                                 sd(PAML_10.0_summary["N0_N1", ]),
                                 sd(PAML_50.0_summary["N0_N1", ]),
                                 sd(PAML_100.0_summary["N0_N1", ]),
                                 sd(PAML_500.0_summary["N0_N1", ]))
PAML.N0.N1.paralog2.mean.list <- c(mean(PAML_3.0_summary["N0_N6", ]), 
                                   mean(PAML_10.0_summary["N0_N6", ]),
                                   mean(PAML_50.0_summary["N0_N6", ]),
                                   mean(PAML_100.0_summary["N0_N6", ]),
                                   mean(PAML_500.0_summary["N0_N6", ]))
PAML.N0.N1.paralog2.sd.list <- c(sd(PAML_3.0_summary["N0_N6", ]), 
                                 sd(PAML_10.0_summary["N0_N6", ]),
                                 sd(PAML_50.0_summary["N0_N6", ]),
                                 sd(PAML_100.0_summary["N0_N6", ]),
                                 sd(PAML_500.0_summary["N0_N6", ]))
# N0_kluyveri
IGC.N0.kluyveri.mean.list <- c(mean(geo_3.0_summary["(N0,kluyveri)", ]), mean(geo_10.0_summary["(N0,kluyveri)", ]),
                               mean(geo_50.0_summary["(N0,kluyveri)", ]), mean(geo_100.0_summary["(N0,kluyveri)", ]),
                               mean(geo_500.0_summary["(N0,kluyveri)", ]))
IGC.N0.kluyveri.sd.list <- c(sd(geo_3.0_summary["(N0,kluyveri)", ]), sd(geo_10.0_summary["(N0,kluyveri)", ]),
                             sd(geo_50.0_summary["(N0,kluyveri)", ]), sd(geo_100.0_summary["(N0,kluyveri)", ]),
                             sd(geo_500.0_summary["(N0,kluyveri)", ]))
PAML.N0.kluyveri.paralog1.mean.list <- c(mean(PAML_3.0_summary["N0_kluyveriYDR418W", ]), 
                                         mean(PAML_10.0_summary["N0_kluyveriYDR418W", ]),
                                         mean(PAML_50.0_summary["N0_kluyveriYDR418W", ]),
                                         mean(PAML_100.0_summary["N0_kluyveriYDR418W", ]),
                                         mean(PAML_500.0_summary["N0_kluyveriYDR418W", ]))
PAML.N0.kluyveri.paralog1.sd.list <- c(sd(PAML_3.0_summary["N0_kluyveriYDR418W", ]), 
                                         sd(PAML_10.0_summary["N0_kluyveriYDR418W", ]),
                                         sd(PAML_50.0_summary["N0_kluyveriYDR418W", ]),
                                         sd(PAML_100.0_summary["N0_kluyveriYDR418W", ]),
                                         sd(PAML_500.0_summary["N0_kluyveriYDR418W", ] ))
# N1_N2
IGC.N1.N2.mean.list <- c(mean(geo_3.0_summary["(N1,N2)", ]), mean(geo_10.0_summary["(N1,N2)", ]),
                         mean(geo_50.0_summary["(N1,N2)", ]), mean(geo_100.0_summary["(N1,N2)", ]),
                         mean(geo_500.0_summary["(N1,N2)", ]))
IGC.N1.N2.sd.list <- c(sd(geo_3.0_summary["(N1,N2)", ]), sd(geo_10.0_summary["(N1,N2)", ]),
                       sd(geo_50.0_summary["(N1,N2)", ]), sd(geo_100.0_summary["(N1,N2)", ]),
                       sd(geo_500.0_summary["(N1,N2)", ]))

PAML.N1.N2.paralog1.mean.list <- c(mean(PAML_3.0_summary["N1_N2", ]), 
                                   mean(PAML_10.0_summary["N1_N2", ]),
                                   mean(PAML_50.0_summary["N1_N2", ]),
                                   mean(PAML_100.0_summary["N1_N2", ]),
                                   mean(PAML_500.0_summary["N1_N2", ]))
PAML.N1.N2.paralog1.sd.list <- c(sd(PAML_3.0_summary["N1_N2", ]), 
                                 sd(PAML_10.0_summary["N1_N2", ]),
                                 sd(PAML_50.0_summary["N1_N2", ]),
                                 sd(PAML_100.0_summary["N1_N2", ]),
                                 sd(PAML_500.0_summary["N1_N2", ]))
PAML.N1.N2.paralog2.mean.list <- c(mean(PAML_3.0_summary["N6_N7", ]), 
                                   mean(PAML_10.0_summary["N6_N7", ]),
                                   mean(PAML_50.0_summary["N6_N7", ]),
                                   mean(PAML_100.0_summary["N6_N7", ]),
                                   mean(PAML_500.0_summary["N6_N7", ]))
PAML.N1.N2.paralog2.sd.list <- c(sd(PAML_3.0_summary["N6_N7", ]), 
                                 sd(PAML_10.0_summary["N6_N7", ]),
                                 sd(PAML_50.0_summary["N6_N7", ]),
                                 sd(PAML_100.0_summary["N6_N7", ]),
                                 sd(PAML_500.0_summary["N6_N7", ]))

# N1_castellii
IGC.N1.castellii.mean.list <- c(mean(geo_3.0_summary["(N1,castellii)", ]), mean(geo_10.0_summary["(N1,castellii)", ]),
                                mean(geo_50.0_summary["(N1,castellii)", ]), mean(geo_100.0_summary["(N1,castellii)", ]),
                                mean(geo_500.0_summary["(N1,castellii)", ]))
IGC.N1.castellii.sd.list <- c(sd(geo_3.0_summary["(N1,castellii)", ]), sd(geo_10.0_summary["(N1,castellii)", ]),
                              sd(geo_50.0_summary["(N1,castellii)", ]), sd(geo_100.0_summary["(N1,castellii)", ]),
                              sd(geo_500.0_summary["(N1,castellii)", ]))

PAML.N1.castellii.paralog1.mean.list <- c(mean(PAML_3.0_summary["N1_castelliiYDR418W", ]), 
                                          mean(PAML_10.0_summary["N1_castelliiYDR418W", ]),
                                          mean(PAML_50.0_summary["N1_castelliiYDR418W", ]),
                                          mean(PAML_100.0_summary["N1_castelliiYDR418W", ]),
                                          mean(PAML_500.0_summary["N1_castelliiYDR418W", ]))
PAML.N1.castellii.paralog1.sd.list <- c(sd(PAML_3.0_summary["N1_castelliiYDR418W", ]), 
                                        sd(PAML_10.0_summary["N1_castelliiYDR418W", ]),
                                        sd(PAML_50.0_summary["N1_castelliiYDR418W", ]),
                                        sd(PAML_100.0_summary["N1_castelliiYDR418W", ]),
                                        sd(PAML_500.0_summary["N1_castelliiYDR418W", ]))

PAML.N1.castellii.paralog2.mean.list <- c(mean(PAML_3.0_summary["N6_castelliiYEL054C", ]), 
                                          mean(PAML_10.0_summary["N6_castelliiYEL054C", ]),
                                          mean(PAML_50.0_summary["N6_castelliiYEL054C", ]),
                                          mean(PAML_100.0_summary["N6_castelliiYEL054C", ]),
                                          mean(PAML_500.0_summary["N6_castelliiYEL054C", ]))
PAML.N1.castellii.paralog2.sd.list <- c(sd(PAML_3.0_summary["N6_castelliiYEL054C", ]), 
                                        sd(PAML_10.0_summary["N6_castelliiYEL054C", ]),
                                        sd(PAML_50.0_summary["N6_castelliiYEL054C", ]),
                                        sd(PAML_100.0_summary["N6_castelliiYEL054C", ]),
                                        sd(PAML_500.0_summary["N6_castelliiYEL054C", ]))
# N2_N3
IGC.N2.N3.mean.list <- c(mean(geo_3.0_summary["(N2,N3)", ]), mean(geo_10.0_summary["(N2,N3)", ]),
                         mean(geo_50.0_summary["(N2,N3)", ]), mean(geo_100.0_summary["(N2,N3)", ]),
                         mean(geo_500.0_summary["(N2,N3)", ]))
IGC.N2.N3.sd.list <- c(sd(geo_3.0_summary["(N2,N3)", ]), sd(geo_10.0_summary["(N2,N3)", ]),
                       sd(geo_50.0_summary["(N2,N3)", ]), sd(geo_100.0_summary["(N2,N3)", ]),
                       sd(geo_500.0_summary["(N2,N3)", ]))

PAML.N2.N3.paralog1.mean.list <- c(mean(PAML_3.0_summary["N2_N3", ]), 
                                   mean(PAML_10.0_summary["N2_N3", ]),
                                   mean(PAML_50.0_summary["N2_N3", ]),
                                   mean(PAML_100.0_summary["N2_N3", ]),
                                   mean(PAML_500.0_summary["N2_N3", ]))
PAML.N2.N3.paralog1.sd.list <- c(sd(PAML_3.0_summary["N2_N3", ]), 
                                 sd(PAML_10.0_summary["N2_N3", ]),
                                 sd(PAML_50.0_summary["N2_N3", ]),
                                 sd(PAML_100.0_summary["N2_N3", ]),
                                 sd(PAML_500.0_summary["N2_N3", ]))
PAML.N2.N3.paralog2.mean.list <- c(mean(PAML_3.0_summary["N7_N8", ]), 
                                   mean(PAML_10.0_summary["N7_N8", ]),
                                   mean(PAML_50.0_summary["N7_N8", ]),
                                   mean(PAML_100.0_summary["N7_N8", ]),
                                   mean(PAML_500.0_summary["N7_N8", ]))
PAML.N2.N3.paralog2.sd.list <- c(sd(PAML_3.0_summary["N7_N8", ]), 
                                 sd(PAML_10.0_summary["N7_N8", ]),
                                 sd(PAML_50.0_summary["N7_N8", ]),
                                 sd(PAML_100.0_summary["N7_N8", ]),
                                 sd(PAML_500.0_summary["N7_N8", ]))
# N2_bayanus
IGC.N2.bayanus.mean.list <- c(mean(geo_3.0_summary["(N2,bayanus)", ]), mean(geo_10.0_summary["(N2,bayanus)", ]),
                              mean(geo_50.0_summary["(N2,bayanus)", ]), mean(geo_100.0_summary["(N2,bayanus)", ]),
                              mean(geo_500.0_summary["(N2,bayanus)", ]))
IGC.N2.bayanus.sd.list <- c(sd(geo_3.0_summary["(N2,bayanus)", ]), sd(geo_10.0_summary["(N2,bayanus)", ]),
                            sd(geo_50.0_summary["(N2,bayanus)", ]), sd(geo_100.0_summary["(N2,bayanus)", ]),
                            sd(geo_500.0_summary["(N2,bayanus)", ]))

PAML.N2.bayanus.paralog1.mean.list <- c(mean(PAML_3.0_summary["N2_bayanusYDR418W", ]), 
                                        mean(PAML_10.0_summary["N2_bayanusYDR418W", ]),
                                        mean(PAML_50.0_summary["N2_bayanusYDR418W", ]),
                                        mean(PAML_100.0_summary["N2_bayanusYDR418W", ]),
                                        mean(PAML_500.0_summary["N2_bayanusYDR418W", ]))
PAML.N2.bayanus.paralog1.sd.list <- c(sd(PAML_3.0_summary["N2_bayanusYDR418W", ]), 
                                      sd(PAML_10.0_summary["N2_bayanusYDR418W", ]),
                                      sd(PAML_50.0_summary["N2_bayanusYDR418W", ]),
                                      sd(PAML_100.0_summary["N2_bayanusYDR418W", ]),
                                      sd(PAML_500.0_summary["N2_bayanusYDR418W", ]))
PAML.N2.bayanus.paralog2.mean.list <- c(mean(PAML_3.0_summary["N7_bayanusYEL054C", ]), 
                                        mean(PAML_10.0_summary["N7_bayanusYEL054C", ]),
                                        mean(PAML_50.0_summary["N7_bayanusYEL054C", ]),
                                        mean(PAML_100.0_summary["N7_bayanusYEL054C", ]),
                                        mean(PAML_500.0_summary["N7_bayanusYEL054C", ]))
PAML.N2.bayanus.paralog2.sd.list <- c(sd(PAML_3.0_summary["N7_bayanusYEL054C", ]), 
                                      sd(PAML_10.0_summary["N7_bayanusYEL054C", ]),
                                      sd(PAML_50.0_summary["N7_bayanusYEL054C", ]),
                                      sd(PAML_100.0_summary["N7_bayanusYEL054C", ]),
                                      sd(PAML_500.0_summary["N7_bayanusYEL054C", ]))

# N3_N4
IGC.N3.N4.mean.list <- c(mean(geo_3.0_summary["(N3,N4)", ]), mean(geo_10.0_summary["(N3,N4)", ]),
                         mean(geo_50.0_summary["(N3,N4)", ]), mean(geo_100.0_summary["(N3,N4)", ]),
                         mean(geo_500.0_summary["(N3,N4)", ]))
IGC.N3.N4.sd.list <- c(sd(geo_3.0_summary["(N3,N4)", ]), sd(geo_10.0_summary["(N3,N4)", ]),
                       sd(geo_50.0_summary["(N3,N4)", ]), sd(geo_100.0_summary["(N3,N4)", ]),
                       sd(geo_500.0_summary["(N3,N4)", ]))
PAML.N3.N4.paralog1.mean.list <- c(mean(PAML_3.0_summary["N3_N4", ]), 
                                   mean(PAML_10.0_summary["N3_N4", ]),
                                   mean(PAML_50.0_summary["N3_N4", ]),
                                   mean(PAML_100.0_summary["N3_N4", ]),
                                   mean(PAML_500.0_summary["N3_N4", ]))
PAML.N3.N4.paralog1.sd.list <- c(sd(PAML_3.0_summary["N3_N4", ]), 
                                 sd(PAML_10.0_summary["N3_N4", ]),
                                 sd(PAML_50.0_summary["N3_N4", ]),
                                 sd(PAML_100.0_summary["N3_N4", ]),
                                 sd(PAML_500.0_summary["N3_N4", ]))
PAML.N3.N4.paralog2.mean.list <- c(mean(PAML_3.0_summary["N8_N9", ]), 
                                   mean(PAML_10.0_summary["N8_N9", ]),
                                   mean(PAML_50.0_summary["N8_N9", ]),
                                   mean(PAML_100.0_summary["N8_N9", ]),
                                   mean(PAML_500.0_summary["N8_N9", ]))
PAML.N3.N4.paralog2.sd.list <- c(sd(PAML_3.0_summary["N8_N9", ]), 
                                 sd(PAML_10.0_summary["N8_N9", ]),
                                 sd(PAML_50.0_summary["N8_N9", ]),
                                 sd(PAML_100.0_summary["N8_N9", ]),
                                 sd(PAML_500.0_summary["N8_N9", ]))
# N3_kudriavzevii
IGC.N3.kudriavzevii.mean.list <- c(mean(geo_3.0_summary["(N3,kudriavzevii)", ]), mean(geo_10.0_summary["(N3,kudriavzevii)", ]),
                                   mean(geo_50.0_summary["(N3,kudriavzevii)", ]), mean(geo_100.0_summary["(N3,kudriavzevii)", ]),
                                   mean(geo_500.0_summary["(N3,kudriavzevii)", ]))
IGC.N3.kudriavzevii.sd.list <- c(sd(geo_3.0_summary["(N3,kudriavzevii)", ]), sd(geo_10.0_summary["(N3,kudriavzevii)", ]),
                                 sd(geo_50.0_summary["(N3,kudriavzevii)", ]), sd(geo_100.0_summary["(N3,kudriavzevii)", ]),
                                 sd(geo_500.0_summary["(N3,kudriavzevii)", ]))
PAML.N3.kudriavzevii.paralog1.mean.list <- c(mean(PAML_3.0_summary["N3_kudriavzeviiYDR418W", ]), 
                                             mean(PAML_10.0_summary["N3_kudriavzeviiYDR418W", ]),
                                             mean(PAML_50.0_summary["N3_kudriavzeviiYDR418W", ]),
                                             mean(PAML_100.0_summary["N3_kudriavzeviiYDR418W", ]),
                                             mean(PAML_500.0_summary["N3_kudriavzeviiYDR418W", ]))
PAML.N3.kudriavzevii.paralog1.sd.list <- c(sd(PAML_3.0_summary["N3_kudriavzeviiYDR418W", ]), 
                                           sd(PAML_10.0_summary["N3_kudriavzeviiYDR418W", ]),
                                           sd(PAML_50.0_summary["N3_kudriavzeviiYDR418W", ]),
                                           sd(PAML_100.0_summary["N3_kudriavzeviiYDR418W", ]),
                                           sd(PAML_500.0_summary["N3_kudriavzeviiYDR418W", ]))
PAML.N3.kudriavzevii.paralog2.mean.list <- c(mean(PAML_3.0_summary["N8_kudriavzeviiYEL054C", ]), 
                                             mean(PAML_10.0_summary["N8_kudriavzeviiYEL054C", ]),
                                             mean(PAML_50.0_summary["N8_kudriavzeviiYEL054C", ]),
                                             mean(PAML_100.0_summary["N8_kudriavzeviiYEL054C", ]),
                                             mean(PAML_500.0_summary["N8_kudriavzeviiYEL054C", ]))
PAML.N3.kudriavzevii.paralog2.sd.list <- c(sd(PAML_3.0_summary["N8_kudriavzeviiYEL054C", ]), 
                                           sd(PAML_10.0_summary["N8_kudriavzeviiYEL054C", ]),
                                           sd(PAML_50.0_summary["N8_kudriavzeviiYEL054C", ]),
                                           sd(PAML_100.0_summary["N8_kudriavzeviiYEL054C", ]),
                                           sd(PAML_500.0_summary["N8_kudriavzeviiYEL054C", ]))
# N4_N5
IGC.N4.N5.mean.list <- c(mean(geo_3.0_summary["(N4,N5)", ]), mean(geo_10.0_summary["(N4,N5)", ]),
                         mean(geo_50.0_summary["(N4,N5)", ]), mean(geo_100.0_summary["(N4,N5)", ]),
                         mean(geo_500.0_summary["(N4,N5)", ]))
IGC.N4.N5.sd.list <- c(sd(geo_3.0_summary["(N4,N5)", ]), sd(geo_10.0_summary["(N4,N5)", ]),
                       sd(geo_50.0_summary["(N4,N5)", ]), sd(geo_100.0_summary["(N4,N5)", ]),
                       sd(geo_500.0_summary["(N4,N5)", ]))
PAML.N4.N5.paralog1.mean.list <- c(mean(PAML_3.0_summary["N4_N5", ]), 
                                   mean(PAML_10.0_summary["N4_N5", ]),
                                   mean(PAML_50.0_summary["N4_N5", ]),
                                   mean(PAML_100.0_summary["N4_N5", ]),
                                   mean(PAML_500.0_summary["N4_N5", ]))
PAML.N4.N5.paralog1.sd.list <- c(sd(PAML_3.0_summary["N4_N5", ]), 
                                 sd(PAML_10.0_summary["N4_N5", ]),
                                 sd(PAML_50.0_summary["N4_N5", ]),
                                 sd(PAML_100.0_summary["N4_N5", ]),
                                 sd(PAML_500.0_summary["N4_N5", ]))
PAML.N4.N5.paralog2.mean.list <- c(mean(PAML_3.0_summary["N9_N10", ]), 
                                   mean(PAML_10.0_summary["N9_N10", ]),
                                   mean(PAML_50.0_summary["N9_N10", ]),
                                   mean(PAML_100.0_summary["N9_N10", ]),
                                   mean(PAML_500.0_summary["N9_N10", ]))
PAML.N4.N5.paralog2.sd.list <- c(sd(PAML_3.0_summary["N9_N10", ]), 
                                 sd(PAML_10.0_summary["N9_N10", ]),
                                 sd(PAML_50.0_summary["N9_N10", ]),
                                 sd(PAML_100.0_summary["N9_N10", ]),
                                 sd(PAML_500.0_summary["N9_N10", ]))
# N4_mikatae
IGC.N4.mikatae.mean.list <- c(mean(geo_3.0_summary["(N4,mikatae)", ]), mean(geo_10.0_summary["(N4,mikatae)", ]),
                              mean(geo_50.0_summary["(N4,mikatae)", ]), mean(geo_100.0_summary["(N4,mikatae)", ]),
                              mean(geo_500.0_summary["(N4,mikatae)", ]))
IGC.N4.mikatae.sd.list <- c(sd(geo_3.0_summary["(N4,mikatae)", ]), sd(geo_10.0_summary["(N4,mikatae)", ]),
                            sd(geo_50.0_summary["(N4,mikatae)", ]), sd(geo_100.0_summary["(N4,mikatae)", ]),
                            sd(geo_500.0_summary["(N4,mikatae)", ]))
PAML.N4.mikatae.paralog1.mean.list <- c(mean(PAML_3.0_summary["N4_mikataeYDR418W", ]), 
                                        mean(PAML_10.0_summary["N4_mikataeYDR418W", ]),
                                        mean(PAML_50.0_summary["N4_mikataeYDR418W", ]),
                                        mean(PAML_100.0_summary["N4_mikataeYDR418W", ]),
                                        mean(PAML_500.0_summary["N4_mikataeYDR418W", ]))
PAML.N4.mikatae.paralog1.sd.list <- c(sd(PAML_3.0_summary["N4_mikataeYDR418W", ]), 
                                      sd(PAML_10.0_summary["N4_mikataeYDR418W", ]),
                                      sd(PAML_50.0_summary["N4_mikataeYDR418W", ]),
                                      sd(PAML_100.0_summary["N4_mikataeYDR418W", ]),
                                      sd(PAML_500.0_summary["N4_mikataeYDR418W", ]))
PAML.N4.mikatae.paralog2.mean.list <- c(mean(PAML_3.0_summary["N9_mikataeYEL054C", ]), 
                                        mean(PAML_10.0_summary["N9_mikataeYEL054C", ]),
                                        mean(PAML_50.0_summary["N9_mikataeYEL054C", ]),
                                        mean(PAML_100.0_summary["N9_mikataeYEL054C", ]),
                                        mean(PAML_500.0_summary["N9_mikataeYEL054C", ]))
PAML.N4.mikatae.paralog2.sd.list <- c(sd(PAML_3.0_summary["N9_mikataeYEL054C", ]), 
                                      sd(PAML_10.0_summary["N9_mikataeYEL054C", ]),
                                      sd(PAML_50.0_summary["N9_mikataeYEL054C", ]),
                                      sd(PAML_100.0_summary["N9_mikataeYEL054C", ]),
                                      sd(PAML_500.0_summary["N9_mikataeYEL054C", ]))
# N5_cerevisiae
IGC.N5.cerevisiae.mean.list <- c(mean(geo_3.0_summary["(N5,cerevisiae)", ]), mean(geo_10.0_summary["(N5,cerevisiae)", ]),
                                 mean(geo_50.0_summary["(N5,cerevisiae)", ]), mean(geo_100.0_summary["(N5,cerevisiae)", ]),
                                 mean(geo_500.0_summary["(N5,cerevisiae)", ]))
IGC.N5.cerevisiae.sd.list <- c(sd(geo_3.0_summary["(N5,cerevisiae)", ]), sd(geo_10.0_summary["(N5,cerevisiae)", ]),
                               sd(geo_50.0_summary["(N5,cerevisiae)", ]), sd(geo_100.0_summary["(N5,cerevisiae)", ]),
                               sd(geo_500.0_summary["(N5,cerevisiae)", ]))
PAML.N5.cerevisiae.paralog1.mean.list <- c(mean(PAML_3.0_summary["N5_cerevisiaeYDR418W", ]), 
                                           mean(PAML_10.0_summary["N5_cerevisiaeYDR418W", ]),
                                           mean(PAML_50.0_summary["N5_cerevisiaeYDR418W", ]),
                                           mean(PAML_100.0_summary["N5_cerevisiaeYDR418W", ]),
                                           mean(PAML_500.0_summary["N5_cerevisiaeYDR418W", ]))
PAML.N5.cerevisiae.paralog1.sd.list <- c(sd(PAML_3.0_summary["N5_cerevisiaeYDR418W", ]), 
                                         sd(PAML_10.0_summary["N5_cerevisiaeYDR418W", ]),
                                         sd(PAML_50.0_summary["N5_cerevisiaeYDR418W", ]),
                                         sd(PAML_100.0_summary["N5_cerevisiaeYDR418W", ]),
                                         sd(PAML_500.0_summary["N5_cerevisiaeYDR418W", ]))
PAML.N5.cerevisiae.paralog2.mean.list <- c(mean(PAML_3.0_summary["N10_cerevisiaeYEL054C", ]), 
                                           mean(PAML_10.0_summary["N10_cerevisiaeYEL054C", ]),
                                           mean(PAML_50.0_summary["N10_cerevisiaeYEL054C", ]),
                                           mean(PAML_100.0_summary["N10_cerevisiaeYEL054C", ]),
                                           mean(PAML_500.0_summary["N10_cerevisiaeYEL054C", ]))
PAML.N5.cerevisiae.paralog2.sd.list <- c(sd(PAML_3.0_summary["N10_cerevisiaeYEL054C", ]), 
                                         sd(PAML_10.0_summary["N10_cerevisiaeYEL054C", ]),
                                         sd(PAML_50.0_summary["N10_cerevisiaeYEL054C", ]),
                                         sd(PAML_100.0_summary["N10_cerevisiaeYEL054C", ]),
                                         sd(PAML_500.0_summary["N10_cerevisiaeYEL054C", ]))
# N5_paradoxus
IGC.N5.paradoxus.mean.list <- c(mean(geo_3.0_summary["(N5,paradoxus)", ]), mean(geo_10.0_summary["(N5,paradoxus)", ]),
                                mean(geo_50.0_summary["(N5,paradoxus)", ]), mean(geo_100.0_summary["(N5,paradoxus)", ]),
                                mean(geo_500.0_summary["(N5,paradoxus)", ]))
IGC.N5.paradoxus.sd.list <- c(sd(geo_3.0_summary["(N5,paradoxus)", ]), sd(geo_10.0_summary["(N5,paradoxus)", ]),
                              sd(geo_50.0_summary["(N5,paradoxus)", ]), sd(geo_100.0_summary["(N5,paradoxus)", ]),
                              sd(geo_500.0_summary["(N5,paradoxus)", ]))
PAML.N5.paradoxus.paralog1.mean.list <- c(mean(PAML_3.0_summary["N5_paradoxusYDR418W", ]), 
                                          mean(PAML_10.0_summary["N5_paradoxusYDR418W", ]),
                                          mean(PAML_50.0_summary["N5_paradoxusYDR418W", ]),
                                          mean(PAML_100.0_summary["N5_paradoxusYDR418W", ]),
                                          mean(PAML_500.0_summary["N5_paradoxusYDR418W", ]))
PAML.N5.paradoxus.paralog1.sd.list <- c(sd(PAML_3.0_summary["N5_paradoxusYDR418W", ]), 
                                        sd(PAML_10.0_summary["N5_paradoxusYDR418W", ]),
                                        sd(PAML_50.0_summary["N5_paradoxusYDR418W", ]),
                                        sd(PAML_100.0_summary["N5_paradoxusYDR418W", ]),
                                        sd(PAML_500.0_summary["N5_paradoxusYDR418W", ]))
PAML.N5.paradoxus.paralog2.mean.list <- c(mean(PAML_3.0_summary["N10_paradoxusYEL054C", ]), 
                                          mean(PAML_10.0_summary["N10_paradoxusYEL054C", ]),
                                          mean(PAML_50.0_summary["N10_paradoxusYEL054C", ]),
                                          mean(PAML_100.0_summary["N10_paradoxusYEL054C", ]),
                                          mean(PAML_500.0_summary["N10_paradoxusYEL054C", ]))
PAML.N5.paradoxus.paralog2.sd.list <- c(sd(PAML_3.0_summary["N10_paradoxusYEL054C", ]), 
                                        sd(PAML_10.0_summary["N10_paradoxusYEL054C", ]),
                                        sd(PAML_50.0_summary["N10_paradoxusYEL054C", ]),
                                        sd(PAML_100.0_summary["N10_paradoxusYEL054C", ]),
                                        sd(PAML_500.0_summary["N10_paradoxusYEL054C", ]))

Now plot all these tree branch length estimates from IGC + MG94 model

# N0_N1
plot(IGC.geo.list, IGC.N0.N1.mean.list)
abline(h = 0.0197240946542)

plot(IGC.geo.list, IGC.N0.N1.sd.list)

# N0_kluyveri
plot(IGC.geo.list, IGC.N0.kluyveri.mean.list)
abline(h = 0.215682181791)

plot(IGC.geo.list, IGC.N0.kluyveri.sd.list)

# N1_N2
plot(IGC.geo.list, IGC.N1.N2.mean.list)
abline(h = 0.20925129872)

plot(IGC.geo.list, IGC.N1.N2.sd.list)

# N1_castellii
plot(IGC.geo.list, IGC.N1.castellii.mean.list)
abline(h = 0.171684721483)

plot(IGC.geo.list, IGC.N1.castellii.sd.list)

# N2_N3
plot(IGC.geo.list, IGC.N2.N3.mean.list)
abline(h = 0.0257112589202)

plot(IGC.geo.list, IGC.N2.N3.sd.list)

# N2_bayanus
plot(IGC.geo.list, IGC.N2.bayanus.mean.list)
abline(h = 0.0266075664688)

plot(IGC.geo.list, IGC.N2.bayanus.sd.list)

# N3_N4
plot(IGC.geo.list, IGC.N3.N4.mean.list)
abline(h = 0.0321083243449)

plot(IGC.geo.list, IGC.N3.N4.sd.list)

# N3_kudriavzevii
plot(IGC.geo.list, IGC.N3.kudriavzevii.mean.list)
abline(h = 0.0853588718458)

plot(IGC.geo.list, IGC.N3.kudriavzevii.sd.list)

# N4_N5
plot(IGC.geo.list, IGC.N4.N5.mean.list)
abline(h = 0.024947887926)

plot(IGC.geo.list, IGC.N4.N5.sd.list)

# N4_mikatae
plot(IGC.geo.list, IGC.N4.mikatae.mean.list)
abline(h = 0.0566627496729)

plot(IGC.geo.list, IGC.N4.mikatae.sd.list)

# N5_cerevisiae
plot(IGC.geo.list, IGC.N5.cerevisiae.mean.list)
abline(h = 0.0581451177847)

plot(IGC.geo.list, IGC.N5.cerevisiae.sd.list)

# N5_paradoxus
plot(IGC.geo.list, IGC.N5.paradoxus.mean.list)
abline(h = 0.0218788166581)

plot(IGC.geo.list, IGC.N5.paradoxus.sd.list)

Feb 1st update

Summarize posterior #IGC + #mutation per site of each branch

post.branch.geo.3.0 <- (geo_3.0_summary[34:45, ] + geo_3.0_summary[46:57, ] + geo_3.0_summary[58:69, ])/geo_3.0_summary[1, ]
row.names(post.branch.geo.3.0) <- row.names(geo_3.0_summary)[10:21]
post.branch.geo.10.0 <- (geo_10.0_summary[34:45, ] + geo_10.0_summary[46:57, ] + geo_10.0_summary[58:69, ])/geo_10.0_summary[1, ]
row.names(post.branch.geo.10.0) <- row.names(geo_10.0_summary)[10:21]
post.branch.geo.50.0 <- (geo_50.0_summary[34:45, ] + geo_50.0_summary[46:57, ] + geo_50.0_summary[58:69, ])/geo_50.0_summary[1, ]
row.names(post.branch.geo.50.0) <- row.names(geo_50.0_summary)[10:21]
post.branch.geo.100.0 <- (geo_100.0_summary[34:45, ] + geo_100.0_summary[46:57, ] + geo_100.0_summary[58:69, ])/geo_100.0_summary[1, ]
row.names(post.branch.geo.100.0) <- row.names(geo_100.0_summary)[10:21]
post.branch.geo.500.0 <- (geo_500.0_summary[34:45, ] + geo_500.0_summary[46:57, ] + geo_500.0_summary[58:69, ])/geo_500.0_summary[1, ]
row.names(post.branch.geo.500.0) <- row.names(geo_500.0_summary)[10:21]
# N0_N1
mean.post.N0.N1.list <- c(mean(post.branch.geo.3.0["(N0,N1)", ]), mean(post.branch.geo.10.0["(N0,N1)", ]),
                          mean(post.branch.geo.50.0["(N0,N1)", ]), mean(post.branch.geo.100.0["(N0,N1)", ]), 
                          mean(post.branch.geo.500.0["(N0,N1)", ])) / 2.0
# N0_kluyveri
mean.post.N0.kluyveri.list <- c(mean(post.branch.geo.3.0["(N0,kluyveri)", ]), mean(post.branch.geo.10.0["(N0,kluyveri)", ]),
                          mean(post.branch.geo.50.0["(N0,kluyveri)", ]), mean(post.branch.geo.100.0["(N0,kluyveri)", ]), 
                          mean(post.branch.geo.500.0["(N0,kluyveri)", ]))

# N1_N2
mean.post.N1.N2.list <- c(mean(post.branch.geo.3.0["(N1,N2)", ]), mean(post.branch.geo.10.0["(N1,N2)", ]),
                          mean(post.branch.geo.50.0["(N1,N2)", ]), mean(post.branch.geo.100.0["(N1,N2)", ]), 
                          mean(post.branch.geo.500.0["(N1,N2)", ])) / 2.0

# N1_castellii
mean.post.N1.castellii.list <- c(mean(post.branch.geo.3.0["(N1,castellii)", ]), mean(post.branch.geo.10.0["(N1,castellii)", ]),
                          mean(post.branch.geo.50.0["(N1,castellii)", ]), mean(post.branch.geo.100.0["(N1,castellii)", ]), 
                          mean(post.branch.geo.500.0["(N1,castellii)", ])) / 2.0

# N2_N3
mean.post.N2.N3.list <- c(mean(post.branch.geo.3.0["(N2,N3)", ]), mean(post.branch.geo.10.0["(N2,N3)", ]),
                          mean(post.branch.geo.50.0["(N2,N3)", ]), mean(post.branch.geo.100.0["(N2,N3)", ]), 
                          mean(post.branch.geo.500.0["(N2,N3)", ]))  / 2.0

# N2_bayanus
mean.post.N2.bayanus.list <- c(mean(post.branch.geo.3.0["(N2,bayanus)", ]), mean(post.branch.geo.10.0["(N2,bayanus)", ]),
                          mean(post.branch.geo.50.0["(N2,bayanus)", ]), mean(post.branch.geo.100.0["(N2,bayanus)", ]), 
                          mean(post.branch.geo.500.0["(N2,bayanus)", ]))  / 2.0

# N3_N4
mean.post.N3.N4.list <- c(mean(post.branch.geo.3.0["(N3,N4)", ]), mean(post.branch.geo.10.0["(N3,N4)", ]),
                          mean(post.branch.geo.50.0["(N3,N4)", ]), mean(post.branch.geo.100.0["(N3,N4)", ]), 
                          mean(post.branch.geo.500.0["(N3,N4)", ])) / 2.0

# N3_kudriavzevii
mean.post.N3.kudriavzevii.list <- c(mean(post.branch.geo.3.0["(N3,kudriavzevii)", ]), mean(post.branch.geo.10.0["(N3,kudriavzevii)", ]),
                          mean(post.branch.geo.50.0["(N3,kudriavzevii)", ]), mean(post.branch.geo.100.0["(N3,kudriavzevii)", ]), 
                          mean(post.branch.geo.500.0["(N3,kudriavzevii)", ])) / 2.0

# N4_N5
mean.post.N4.N5.list <- c(mean(post.branch.geo.3.0["(N4,N5)", ]), mean(post.branch.geo.10.0["(N4,N5)", ]),
                          mean(post.branch.geo.50.0["(N4,N5)", ]), mean(post.branch.geo.100.0["(N4,N5)", ]), 
                          mean(post.branch.geo.500.0["(N4,N5)", ])) / 2.0

# N4_mikatae
mean.post.N4.mikatae.list <- c(mean(post.branch.geo.3.0["(N4,mikatae)", ]), mean(post.branch.geo.10.0["(N4,mikatae)", ]),
                          mean(post.branch.geo.50.0["(N4,mikatae)", ]), mean(post.branch.geo.100.0["(N4,mikatae)", ]), 
                          mean(post.branch.geo.500.0["(N4,mikatae)", ])) / 2.0

# N5_cerevisiae
mean.post.N5.cerevisiae.list <- c(mean(post.branch.geo.3.0["(N5,cerevisiae)", ]), mean(post.branch.geo.10.0["(N5,cerevisiae)", ]),
                          mean(post.branch.geo.50.0["(N5,cerevisiae)", ]), mean(post.branch.geo.100.0["(N5,cerevisiae)", ]), 
                          mean(post.branch.geo.500.0["(N5,cerevisiae)", ])) / 2.0

# N5_paradoxus
mean.post.N5.paradoxus.list <- c(mean(post.branch.geo.3.0["(N5,paradoxus)", ]), mean(post.branch.geo.10.0["(N5,paradoxus)", ]),
                          mean(post.branch.geo.50.0["(N5,paradoxus)", ]), mean(post.branch.geo.100.0["(N5,paradoxus)", ]), 
                          mean(post.branch.geo.500.0["(N5,paradoxus)", ])) / 2.0
# N0_N1
sd.post.N0.N1.list <- c(sd(post.branch.geo.3.0["(N0,N1)", ]), sd(post.branch.geo.10.0["(N0,N1)", ]),
                          sd(post.branch.geo.50.0["(N0,N1)", ]), sd(post.branch.geo.100.0["(N0,N1)", ]), 
                          sd(post.branch.geo.500.0["(N0,N1)", ])) / 2.0
# N0_kluyveri
sd.post.N0.kluyveri.list <- c(sd(post.branch.geo.3.0["(N0,kluyveri)", ]), sd(post.branch.geo.10.0["(N0,kluyveri)", ]),
                          sd(post.branch.geo.50.0["(N0,kluyveri)", ]), sd(post.branch.geo.100.0["(N0,kluyveri)", ]), 
                          sd(post.branch.geo.500.0["(N0,kluyveri)", ]))

# N1_N2
sd.post.N1.N2.list <- c(sd(post.branch.geo.3.0["(N1,N2)", ]), sd(post.branch.geo.10.0["(N1,N2)", ]),
                          sd(post.branch.geo.50.0["(N1,N2)", ]), sd(post.branch.geo.100.0["(N1,N2)", ]), 
                          sd(post.branch.geo.500.0["(N1,N2)", ])) / 2.0

# N1_castellii
sd.post.N1.castellii.list <- c(sd(post.branch.geo.3.0["(N1,castellii)", ]), sd(post.branch.geo.10.0["(N1,castellii)", ]),
                          sd(post.branch.geo.50.0["(N1,castellii)", ]), sd(post.branch.geo.100.0["(N1,castellii)", ]), 
                          sd(post.branch.geo.500.0["(N1,castellii)", ])) / 2.0

# N2_N3
sd.post.N2.N3.list <- c(sd(post.branch.geo.3.0["(N2,N3)", ]), sd(post.branch.geo.10.0["(N2,N3)", ]),
                          sd(post.branch.geo.50.0["(N2,N3)", ]), sd(post.branch.geo.100.0["(N2,N3)", ]), 
                          sd(post.branch.geo.500.0["(N2,N3)", ])) / 2.0

# N2_bayanus
sd.post.N2.bayanus.list <- c(sd(post.branch.geo.3.0["(N2,bayanus)", ]), sd(post.branch.geo.10.0["(N2,bayanus)", ]),
                          sd(post.branch.geo.50.0["(N2,bayanus)", ]), sd(post.branch.geo.100.0["(N2,bayanus)", ]), 
                          sd(post.branch.geo.500.0["(N2,bayanus)", ])) / 2.0

# N3_N4
sd.post.N3.N4.list <- c(sd(post.branch.geo.3.0["(N3,N4)", ]), sd(post.branch.geo.10.0["(N3,N4)", ]),
                          sd(post.branch.geo.50.0["(N3,N4)", ]), sd(post.branch.geo.100.0["(N3,N4)", ]), 
                          sd(post.branch.geo.500.0["(N3,N4)", ])) / 2.0

# N3_kudriavzevii
sd.post.N3.kudriavzevii.list <- c(sd(post.branch.geo.3.0["(N3,kudriavzevii)", ]), sd(post.branch.geo.10.0["(N3,kudriavzevii)", ]),
                          sd(post.branch.geo.50.0["(N3,kudriavzevii)", ]), sd(post.branch.geo.100.0["(N3,kudriavzevii)", ]), 
                          sd(post.branch.geo.500.0["(N3,kudriavzevii)", ])) / 2.0

# N4_N5
sd.post.N4.N5.list <- c(sd(post.branch.geo.3.0["(N4,N5)", ]), sd(post.branch.geo.10.0["(N4,N5)", ]),
                          sd(post.branch.geo.50.0["(N4,N5)", ]), sd(post.branch.geo.100.0["(N4,N5)", ]), 
                          sd(post.branch.geo.500.0["(N4,N5)", ])) / 2.0

# N4_mikatae
sd.post.N4.mikatae.list <- c(sd(post.branch.geo.3.0["(N4,mikatae)", ]), sd(post.branch.geo.10.0["(N4,mikatae)", ]),
                          sd(post.branch.geo.50.0["(N4,mikatae)", ]), sd(post.branch.geo.100.0["(N4,mikatae)", ]), 
                          sd(post.branch.geo.500.0["(N4,mikatae)", ])) / 2.0

# N5_cerevisiae
sd.post.N5.cerevisiae.list <- c(sd(post.branch.geo.3.0["(N5,cerevisiae)", ]), sd(post.branch.geo.10.0["(N5,cerevisiae)", ]),
                          sd(post.branch.geo.50.0["(N5,cerevisiae)", ]), sd(post.branch.geo.100.0["(N5,cerevisiae)", ]), 
                          sd(post.branch.geo.500.0["(N5,cerevisiae)", ])) / 2.0

# N5_paradoxus
sd.post.N5.paradoxus.list <- c(sd(post.branch.geo.3.0["(N5,paradoxus)", ]), sd(post.branch.geo.10.0["(N5,paradoxus)", ]),
                          sd(post.branch.geo.50.0["(N5,paradoxus)", ]), sd(post.branch.geo.100.0["(N5,paradoxus)", ]), 
                          sd(post.branch.geo.500.0["(N5,paradoxus)", ])) / 2.0

Now plot same branch length estimates from each paralog in paml results and posterior branch length from IGC + MG94 model

# N0_N1
matplot(IGC.geo.list, cbind(PAML.N0.N1.paralog1.mean.list, PAML.N0.N1.paralog2.mean.list, mean.post.N0.N1.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "mean N0.N1 estimate")
abline(h = 0.0197240946542)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

matplot(IGC.geo.list, cbind(PAML.N0.N1.paralog1.sd.list, PAML.N0.N1.paralog2.sd.list, sd.post.N0.N1.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "sd N0.N1 estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

# N0_kluyveri
matplot(IGC.geo.list, cbind(PAML.N0.kluyveri.paralog1.mean.list, mean.post.N0.kluyveri.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "mean N0.kluyveri estimate")
abline(h = 0.215682181791)
legend("right", legend = c("paralog1", "posterior IGC"), col = 1:3, pch = 1)

matplot(IGC.geo.list, cbind(PAML.N0.kluyveri.paralog1.sd.list, sd.post.N0.kluyveri.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "sd N0.kluyveri estimate")
legend("right", legend = c("paralog1", "posterior IGC"), col = 1:3, pch = 1)

# N1_N2
matplot(IGC.geo.list, cbind(PAML.N1.N2.paralog1.mean.list, PAML.N1.N2.paralog2.mean.list, mean.post.N1.N2.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "mean N1.N2 estimate")
abline(h = 0.20925129872)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

matplot(IGC.geo.list, cbind(PAML.N1.N2.paralog1.sd.list, PAML.N1.N2.paralog2.sd.list, sd.post.N1.N2.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "sd N1.N2 estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

# N1_castellii
matplot(IGC.geo.list, cbind(PAML.N1.castellii.paralog1.mean.list, PAML.N1.castellii.paralog2.mean.list, mean.post.N1.castellii.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "mean N1.castellii estimate")
abline(h = 0.171684721483)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

matplot(IGC.geo.list, cbind(PAML.N1.castellii.paralog1.sd.list, PAML.N1.castellii.paralog2.sd.list, sd.post.N1.castellii.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "sd N1.castellii estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

# N2_N3
matplot(IGC.geo.list, cbind(PAML.N2.N3.paralog1.mean.list, PAML.N2.N3.paralog2.mean.list, mean.post.N2.N3.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "mean N2.N3 estimate")
abline(h = 0.0257112589202)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

matplot(IGC.geo.list, cbind(PAML.N2.N3.paralog1.sd.list, PAML.N2.N3.paralog2.sd.list, sd.post.N2.N3.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "sd N2.N3 estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

# N2_bayanus
matplot(IGC.geo.list, cbind(PAML.N2.bayanus.paralog1.mean.list, PAML.N2.bayanus.paralog2.mean.list, mean.post.N2.bayanus.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "mean N2.bayanus estimate")
abline(h = 0.0266075664688)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

matplot(IGC.geo.list, cbind(PAML.N2.bayanus.paralog1.sd.list, PAML.N2.bayanus.paralog2.sd.list, sd.post.N2.bayanus.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "sd N2.bayanus estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

# N3_N4
matplot(IGC.geo.list, cbind(PAML.N3.N4.paralog1.mean.list, PAML.N3.N4.paralog2.mean.list, mean.post.N3.N4.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "mean N3.N4 estimate")
abline(h = 0.0321083243449)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

matplot(IGC.geo.list, cbind(PAML.N3.N4.paralog1.sd.list, PAML.N3.N4.paralog2.sd.list, sd.post.N3.N4.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "sd N3.N4 estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

# N3_kudriavzevii
matplot(IGC.geo.list, cbind(PAML.N3.kudriavzevii.paralog1.mean.list, PAML.N3.kudriavzevii.paralog2.mean.list, mean.post.N3.kudriavzevii.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "mean N3.kudriavzevii estimate")
abline(h = 0.0853588718458)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

matplot(IGC.geo.list, cbind(PAML.N3.kudriavzevii.paralog1.sd.list, PAML.N3.kudriavzevii.paralog2.sd.list, sd.post.N3.kudriavzevii.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "sd N3.kudriavzevii estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

# N4_N5
matplot(IGC.geo.list, cbind(PAML.N4.N5.paralog1.mean.list, PAML.N4.N5.paralog2.mean.list, mean.post.N4.N5.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "mean N4.N5 estimate")
abline(h = 0.024947887926)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

matplot(IGC.geo.list, cbind(PAML.N4.N5.paralog1.sd.list, PAML.N4.N5.paralog2.sd.list, sd.post.N4.N5.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "sd N4.N5 estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

# N4_mikatae
matplot(IGC.geo.list, cbind(PAML.N4.mikatae.paralog1.mean.list, PAML.N4.mikatae.paralog2.mean.list, mean.post.N4.mikatae.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "mean N4.mikatae estimate")
abline(h = 0.0566627496729)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

matplot(IGC.geo.list, cbind(PAML.N4.mikatae.paralog1.sd.list, PAML.N4.mikatae.paralog2.sd.list, sd.post.N4.mikatae.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "sd N4.mikatae estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

# N5_cerevisiae
matplot(IGC.geo.list, cbind(PAML.N5.cerevisiae.paralog1.mean.list, PAML.N5.cerevisiae.paralog2.mean.list, mean.post.N5.cerevisiae.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "mean N5.cerevisiae estimate")
abline(h = 0.0581451177847)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

matplot(IGC.geo.list, cbind(PAML.N5.cerevisiae.paralog1.sd.list, PAML.N5.cerevisiae.paralog2.sd.list, sd.post.N5.cerevisiae.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "sd N5.cerevisiae estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

# N5_paradoxus
matplot(IGC.geo.list, cbind(PAML.N5.paradoxus.paralog1.mean.list, PAML.N5.paradoxus.paralog2.mean.list, mean.post.N5.paradoxus.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "mean N5.paradoxus estimate")
abline(h = 0.0218788166581)
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

matplot(IGC.geo.list, cbind(PAML.N5.paradoxus.paralog1.sd.list, PAML.N5.paradoxus.paralog2.sd.list, sd.post.N5.paradoxus.list),
        type = c("p"), pch = 1, col = 1:3, ylab = "sd N5.paradoxus estimate")
legend("right", legend = c("paralog1", "paralog2", "posterior IGC"), col = 1:3, pch = 1)

Now save workspace.

#save.image("/Users/xji3/GitFolders/IGCCodonSimulation/All.RData")